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Abstract. Finite volume schemes are commonly used to construct approximate solu- 
tions to conservation laws. In this study we extend the framework of the finite volume 
methods to dispersive water wave models, in particular to Boussinesq type systems. We 
focus mainly on the application of the method to bidirectional nonlinear, dispersive wave 
propagation in one space dimension. Special emphasis is given to important nonlinear 
phenomena such as solitary waves interactions, dispersive shock wave formation and the 
runup of breaking and non-breaking long waves. 



1. Introduction 

The simulation of water waves in realistic and complex environments is a very challenging 
problem. Most of the applications arise from the areas of coastal and naval engineering, 
but also from natural hazards assessment. These applications may require the computation 
of the wave generation [DD07, KDD07], propagation [TG97], interaction with solid bodies, 
the computation of long wave runup [TS94, TS98] and even the extraction of the wave 
energy [SimSl]. Issues like wave breaking, robustness of the numerical algorithm in wet- 
dry processes along with the validity of the mathematical models in the near-shore zone 
are some basic problems in this direction [HP79]. During past years, the classical shallow 
water equations have been employed to solve some of these problems [AC99, DPDIO]: 

{Hu)t + [Hu' + ^H^)^ = gHD,, ^^-^^ 

where H{x, t) := ri{x, t) + D{x) is the total water depth, D{x) describes the depth below the 
mean sea level while ri{x, t) is the free surface elevation, t) denotes the depth-averaged 
fluid velocity and g is the gravity acceleration constant. Mathematically, equations (1.1) 
represent a system of conservation laws describing the propagation of infinitely long waves 
with a hydrostatic pressure assumption. The wave breaking phenomenon is commonly 
assimilated to the formation of shock waves (or hydraulic jumps) which is a common 
feature of hyperbolic p.d.e's. Consequently, the finite volume(FV) method has become 
the method of choice for these problems due to its excellent intrinsic conservative and 
shock-capturing properties [AC99, DK03, DKK08, DPDIO]. Furthermore the shallow water 
equations have been proven in practice to predict accurately the maximum runup of long 
waves [HH70, Syn87, TS94, TS98, KS98, CFVCP04]. 
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On the other hand, various studies have shown that the inclusion of dispersive effects 
is beneficial for the description of long wave propagation and runup processes [Zel91, 
WKGS95, LR02]. Moreover, J. A. Zelt [Zel91] reported a divergence in the prediction 
of the rundown and in the prediction of the reflected wave-train after the wave climbing on 
the shore when a dispersionless model is employed. According to J. A. Zelt, the results of 
the nonlinear dispersive model considered in [Zel91] showed better performance compared 
to (1.1). During the last fifty years numerous dispersive models have been proposed for 
the simulation of long waves [Ser53, Per67, BS76, Nwo93, KCKDOO, MBS03, Mit09]. 

In this work we will study numerically bidirectional water wave models. Specifically, we 
consider the following family of Boussinesq type systems of water wave theory, introduced 
in [BCS02], written in nondimensional, unsealed variables 

+ + M, + = 0, 

Ut + r]x + uux + cq^xx - du^xt = 0, 

where a, 6, c, G M, = ri{x, t), u = u{x, t) are real functions defined for x G M and t > 0. 

For more realistic situations we introduce a new Boussinesq type system with variable 
bottom topography based on Peregrine's system, [Per67]. The new system incorporates a 
very important property — the invariance under vertical translations, thus more appropri- 
ate for practical applications such as wave runup on non-uniform shores. In dimensional 
variables the model reads 



Ht + Qx = 0, 

Qt + {^ + IH% - \^Q..t - {\Hl - \HH^,)Qt + \HH^qJ = gHD,, ^^'^^ 



where H{x, t) = rj{x, t) + D{x), Q{x, t) = H{x, t)u{x, t). 

There is a wide range of numerical methods in the literature for computing approximate 
solutions to these models. Finite difference (FD) schemes [KWC+98, JFBM06, HMK+09], 
finite element methods [BDM07, Mit09, DM08, AVSS09] and spectral methods [PDOl, 
NguOS] have been proposed. More contemporary discontinuous Galerkin (DG) schemes 
have also been adapted with some success to dispersive wave equations [YS02, LSY04, 
ES06, ESB06] while the application of Finite Volumes (FV) or hybrid FV/FD methods 
remain most infrequent for this type of problems. To our knowledge, only a few very recent 
works are in this direction [BBOl, EIK05, BS08, TP09, SM09]. 

Finite volume method is well known for its accuracy, efficiency and robustness for ap- 
proximating solutions to conservation laws and in particular to nonlinear shallow water 
equations (1.1). The aforementioned bidirectional models (1.2) and (1.3) are rewritten in 
a conservative form and discretization by the finite volume method follows. Three different 
numerical fluxes are employed 

• a simple average flux (m-scheme), 

• a central flux, (KT-scheme) [NT90, KTOO], as a representative of central schemes, 

• a characteristic flux (CF-scheme), as a representative of the linearized Riemann 
solvers, [GKC96], 



FINITE VOLUME SCHEMES FOR DISPERSIVE WAVES 



3 



along with TVD, UNO and WENO reconstruction techniques, [Swe84, H087, LOC94]. 
Time discretization is based on Runge-Kutta (RK) methods which preserve the total vari- 
ation diminishing(TVD) property of the finite volume scheme, [S088, GSTOl, SR02]. We 
use explicit RK methods since we work with BBM type systems (1.2) and not with KdV- 
type systems which is well known to be notoriously stiff, [Mit09]. These methods have been 
studied thoroughly in the case of nonlinear conservation laws. The average flux although 
is known to be unstable for conservation laws is proved to be very accurate for nonlinear 
dispersive waves. On the other hand flnite volume methods based on the central flux as 
well as on characteristic flux work equally well for the numerical simulation of waves even 
in realistic environments. 

The performance of the flnite volume method applied to models (1.2) and to the new 
system (1.3) is studied in a systematic way through a series of numerical experiments. Our 
main focus is the evaluation, in terms of accuracy, efficiency and robustness of second order 
finite volume methods compared to high order schemes. In particular, in this study we 
take up on the following points 

• accuracy of the finite volume method in the propagation of solitary waves with very 
satisfactory results. 

• conservation of various invariant quantities during the formation of dispersive shocks 
is studied numerically. The finite element as well as spectral methods break down 
for these experiments. The finite volume method provides very accurate results. 

• interactions of solitary waves are computed with high accuracy. It is shown numer- 
ically that Boussinesq type systems describe better overtaking collisions of solitary 
waves than unidirectional models like KdV-BBM. We compare our results, when- 
ever possible, with experimental measurements with very good agreement. 

• finite volume method allows to use appropriate techniques to treat the transition 
from wet to dry regions and vice versa. These techniques are applied successfully to 
systems with dispersive terms modeling runup of long waves. On the other hand, 
when the model fails due to wave breaking, the method allows to use locally the 
nonlinear shallow water system, thus enabling us to resolve a wide spectrum of 
hydrodynamic phenomena using a single computational framework. 

• it is shown numerically the advantage of using dispersive models over standard 
nonlinear shallow water equations in computing the wave runup and, in particular, 
in capturing the refiected wave. It's also illustrated by an example the importance 
of the system being invariant under vertical translations. 

The paper is organized as follows. In Section 2 Boussinesq type systems are presented 
along with some of their basic properties. A new system with uneven bottom and invariant 
under vertical translations is derived. In Section 3 the finite volume method is presented 
for a general framework incorporating all models. 

Section 4 presents a series of numerical experiments for the Boussinesq systems (1.2). In 
this mathematical setting we validate the finite volume method and measure its accuracy. 
We study the propagation as well as the interaction of solitary waves: we consider in par- 
ticular head-on and overtaking collisions, but also we present results concerning the small 
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dispersion effect. Finally, in Section 5 the new system with variable bottom, (1.3) is stud- 
ied. Numerical simulations of non-breaking and breaking long wave runup are presented 
and compared with experimental data. 

2. Mathematical models 

We present briefly the mathematical models being considered and some of their main 
properties. 

2.1. Dispersive models with flat bottom. We consider the following family of Boussi- 
nesq type systems of water wave theory, introduced in [BCS02], which may be written in 
nondimensional, unsealed variables 

Vt + «x + {vu)x + au^xx - br]^xt = 0, ^2 i) 

Ut + r]x + uu^ + cq^n^n, - du^^t = 0, 

where t] = ri{x,t), u = u{x,t) are real functions defined for x G M and t > 0. Coefficients 
a, b, c and d are defined as 

a = ^(^' - b = \{e' v\ c = 1(1 - e^)y., ^ = ^(1 " " /^)' (2-2) 

where < ^ < 1 and yU, z/ G M. The variables in (2.1) are non-dimensional and unsealed: 
X and t are proportional to position along the channel and time, respectively, while rjlxjt) 
and u{x, t) are proportional to the deviation of the free surface above an undisturbed level 
and to the horizontal velocity of the fluid at a height y = —l + 9{l + rj{x,t)), respectively. 
In terms of these variables the channel bottom is located at y = — 1, (^ = 0), while the free 
surface corresponds to 6* = 1. Boussinesq systems (2.1) with b = d conserve the energy 
functional: 

^i(^) = {v'^{x,t) + {l + ri{x,t))u'^{x,t) - cril{x,t) - aul{x,t)) dx, (2.3) 
Jr 

i.e. Ii{t) = /i(0) for t >0. System (2.1) is derived under the following assumptions: 

e ■= A/ho < 1 , a:= h^/X < 1 S := = 0(1) , 

where 5" is the Stokes (or Ursell) number, A is a typical wave amplitude of fluid of depth ho 
and A is a characteristic wavelength. If one takes S = 1 and switches to scaled, dimension- 
less variables, one may derive from Euler equations a scaled version of (2.1) by appropriate 
asymptotic expansion in powers of e, cf. [BCS02]: 

rjt + + e{riu)^ + e[au^^^ - brj^^t] = ^(e^), 4) 
ut + Vx + euu^ + elcq^xx - du^xt] =0{e^), 

from which we obtain (2.1) by unsealing and neglecting higher order terms 0{e'^). 

We list several examples of particular Boussinesq systems of the form (2.1) that we will 
refer to in the sequel. The initial-value problem for all these systems has been shown to 
be at least nonlinearly well-posed locally in time, cf. [BCS04]. 
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(i) The 'classical' Boussinesq system (/i = 0, 6"^ = 1/3, i.e. a = b = c = 0, d = 1/3 in 
(2.1)), whose initial-value problem is globally well-posed, [Ami84, SchSl], 

Vt + + ivu)x = 0, ^2 5) 

ut + Vx + uux - luxxt = 0. 

(ii) The BBM-BBM system = ^ = 0, 9^ = 2/3, i.e. a = c = 0, b = d = 1/6 in (2.1)), 
whose initial- value problem is locally well-posed, [BC98], 

Vt + Ux + {vu)x - IVxxt = 0, ^2 5) 

Ut + r]x + uUx - \uxxt = 0. 

(iii) The Bona-Smith system (z/ = 0, /i = (4 - Qe'^)/?,{1 - 6^), i.e. a = 0, b = d = 
(3^2 _ i)/5^ c = (2 - 36^) /3, 2/3 < < 1 in (2.1)), whose initial-value problem 
is globally well-posed, cf. [BS76]. The limiting form of this system as 6' — )■ 1, 
corresponding to a = 0, b = d = 1/3, c = —1/3, is the system actually studied by 
Bona and Smith, [BS76]. These systems are given by 

'nt + Ux + {rju)^ - ^^Y^rj^xt = 0, ^2 7) 

Ut + r]x + UUx + ^^^Vxxx - ^^-^Uxxt = 0. 

The existence of solitary wave solutions to the above systems, in some cases the uniqueness 
also, has been proved in [CheOO, Che98, DM04] and in the case of the Bona-Smith type 
systems (2.7), for each 9^ G (7/9, 1), there exists one solitary wave in closed form, [Che98] 

^(0 = ^osech^(AO, /..ON 
u{0 = BviO, ^ ' 



with 

_ 9 6>^-7/9 _ 4(0^-2/3) 

^° ~ 2 ■ 1-62 , Cs - ^2ii -e^)ie^_y Y)-' 



2Y (02-l/3)(e2-2/3) ' V 62-1/3- 



(2.9) 



2.2. Dispersive models with variable bottom. For more realistic applications one 
should consider Boussinesq systems with variable bottom. In his pioneering work Peregrine, 
[Per67], derived the following Boussinesq type system 

„. + [(fl + „H. = o, 

Ut + grix + uu^ - ^[Du)^^t + —Uxxt = 0, 

where ri{x,t) and u{x,t) are defined as before and D{x) describes the water depth below 
its rest position. Many other systems have been derived also, including systems with 
improved dispersion characteristics [Nwo93], high-order Boussinesq systems [MBS03] and 
other generalizations of (2.1), cf. [Mit09]. Most of these systems break Galilean invariance 
and the invariance under vertical translations. This is a restrictive property especially in 
the studies of realistic problems like the water wave runup on non-uniform shores. We note 
also that the complete water wave problem possesses these symmetries, [B082]. 
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To overcome this deficiency we develop a new system, analogous to the original Pere- 
grine's system, [Per67], which is invariant under vertical translations. To derive the system 
we begin with (2.10) written in dimensionless scaled variables (in analogy with (2.4)) 

r]t + [{D + e7])ul = 0, 



Ut + Vx + £UUx - cr' 



2 i,-^^) XXt R 



U. 



xxt 



0(e\ea^) 



(2.11) 



Then by setting H = D + erj, we obtain 

Ht + e{Hu), = 0, 
{Hu)t + {eHu" + 



-a 



HP 
2 



■ Du 
. D 



' xxt 



(2.12) 



Observing that 
we have that 



^)^^ = [2^-^^](Du)-2^^{Du)^ + ^{Du)^^ and that H = D + 0{e) 



Ht + e{Hu)^ = 0, 
{Hu)t + {eHu^ + j-H 



—a 



= \HD, + 0{e\ea^). 

By setting Q = Hu, and using again the relation H 

Ht + eQ^ = 0, 
Qt + {eQ^/H 



(2.13) 



D + 0{e) we have 



-cr 



xxt 



¥1 



lHHxx)Qt + \HHxQxt = \HDj, + 0{e^, ea 



(2.14) 



In dimensional variables, neglecting the higher order terms at the right-hand side we obtain 
Ht + Qx = 0, 

Qt + {Q'/H + IH% - P{H, Q) = gHD,. (2.15) 
P{H, Q) = ^Qxxt — {\HI — \HHxx)Qt + \HHxQxt 

where H{x,t) = ri{x,t) + D{x), Q{x,t) = H{x,t)u{x,t). We underline that system (2.15) 
is invariant under vertical translations and therefore more appropriate for studying long 
wave runup. Moreover, the linearization of the system (2.15) coincides with the original 
Peregrine's system (2.11) and therefore, inherits all its linear dispersive characteristics. On 
the other hand system (2.15) cannot be regarded as a correct asymptotic model to the Euler 
equations since it contains terms of the order 0{ea^). However, such terms considered in 
the correct (small amplitude and long wave) regime are negligible, hence their contribution 
will be negligible. Finally we note that ignoring the dispersive terms P{H, Q) of system 
(2.15) we obtain the shallow water equations (1.1). 

We also note that even though Boussinesq systems are not valid in the near-shore re- 
gion, in practice they appear to predict well the behavior of small amplitude waves from 
moderately deep to shallower waters and for smooth flows, cf. [Zel91]. Of course, more 
accurate systems in the near-shore zone have been derived such as the Serre equations 
(sometimes referred also as Green-Naghdi equations), cf. [Ser53, LB09, DM10]. These 
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systems appeared in practice to model better the breaking phenomena in the near shore 
zone but recent numerical studies of the Serre system showed that unphysical oscillations 
might appear in analogy with the Boussinesq equations during the wave breaking and the 
runup process, [CBB06, CBB07]. 

2.3. Source terms. Nonlinear shallow water model (1.1) and Boussinesq system (1.3) 
may be completed to take into account some dissipative or friction effects which are very 
beneficial in describing the wave breaking phenomena. Usually this is accomplished by 
including appropriate source or dissipative terms into momentum conservation equations 
(1.1) or (1.3). Possible choices are the following : 

u\u\ 

Friction: F{u, H) = -Cm 9 j^ys^ (2-^^) 

Viscosity: V{u, H) = /i— \H—\ , (2.17) 

where is the Manning roughness coefficient and yU denotes the kinematic viscosity of 
the fluid. The particular form of the source terms is suggested by empirical laws, which 
are generally obtained for steady state flows. Similar models have been derived from then 
Navier-Stokes system for incompressible flows with a free surface. More complex friction 
laws can be also formulated to model bottom rugosity effects, etc. 



3. Numerical schemes 

In the present section we generalize the finite volume method to systems (1.2) and (1.3) 
of dispersive PDEs. In our work we rely on corresponding schemes for conservation laws. 
Next we present briefly the finite volume framework for conservation laws. Based on this 
framework we introduce finite volume schemes for the dispersive models. 

3.1. Finite volume method for conservation laws. We consider the initial value prob- 
lem 

wt + Fiw)^ = S(w), X ER,t>0 
w[x,^) = Wo{x), 

where w{x,t) is the state variable, F denotes the fiux and S is the source term. Let 
7" = {xi}, i E Z denotes a partition of M into cells Ci = where = (a^j+i + 

Xi_i)/2 denotes the midpoint of Cj. Let Axi = x^^i — x^_i be the length of the cell 
Cj, Ax^j^i = Xj+i — Xj. Without loss of generality we assume a uniform partition T, 
that is Axi = Ax^^i = Ax, i E Z. Let Wi denotes the cell average of w on Ci i.e 
Wi{t) = w{x,t) dx. Then a simple integration of (3.1) over a cell Ci yields 

;W,{t) + ^[F{w{x,^.,t))-F{w{x,_._,t)))= — 1^ S{w{x,t))dx. (3.2) 



dt 
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3.1.1. Semidis Crete schemes. We now define the semidiscrete finite volume approximation 
of w{x,t). Let xCi denotes tlie cliaracteristic function of tlie cell Ci, we seek a piecewise 
constant function Wh{x,t) = ^jg^ VFj(t)xc,(a;) with 

where J^i+i = -^(^^j+i > W^j+i) is an approximation to , t)) while Si approximates 

the source term Si = Si{Wi) ~ 4- f„ S{w{x,t)) dx. The values W^i,W^i are approxi- 
mations to the point value , t) from cells Cj, Cj+i respectively and J-" is a numerical 

flux function which is consistent and monotone. The values are computed by 

a reconstruction process described below (see Section 3.1.3). 

3.1.2. The numerical fluxes. There are many possible choices for the numerical flux func- 
tion J-". In the present study we choose to work with three following fluxes 

J-(PF, V) = F( , (3.4) 



^^^(H/, r) = 1 {[FiV) + F{W)] - AiW, V)[V-W]}, (3.5) 

:F^P(W, V) = ]^ {[F{V) + F{W)] - A{W, V) [F{V) - F{W)]} . (3.6) 

The average flux (3.4) is the simplest one. It is well known that although this flux is 
unstable for nonlinear conservation laws, it is proven very stable and accurate for nonlinear 
dispersive models. 

The central flux (3.5) is a Lax-Friedrichs type flux and is a representative of central 
schemes [KTOO, NT90]. The operator A is related to the characteristic speeds of the flow 
and is defined as 

A{W, V) = max [p {DF{W)) , p {DF{V))] , (3.7) 

where DF denotes the Jacobian matrix and p{A) is the spectral radius of A. 

The characteristic flux function (3.6), [GKC96, GKCOl], is similar to the upwind flux 
and the operator A in this case is defined by 

A{W, V) = sign (^DF (^^^) ) . (3.8) 

3.1.3. The reconstruction process. The values W^i, are approximations to w{x^_^i,t) 

from cells Ci and Cj+i respectively. The simplest possible choice is to take the piecewise 
constant approximation in each cell, 

W,^, = W,, W^1^ = iy,+i. (3.9) 
The resulting semidiscrete finite volume scheme is formally first order accurate in space. 
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To construct a higher order scheme in space, the piecewise constant data is replaced by 
a piecewise polynomial representation. The main idea here is to construct higher order 
approximations to w{Xi_^_i,t) using the computed cell averages Wi. For this purpose the 
classical MUSCL type (TVD2) linear reconstruction [Kol75, vL79] as well as UN02, [H087] 
or WENO type reconstructions, [LOC94], have been developed. 

The classical TVD2 linear reconstruction is given by the following formulas: 

^^^1 =W, + ^0(r,)(H^,+i - Wi), W^1^ = W,+, - ^0(r,+i)(H^,+2 - W^m), (3.10) 

where rj = , and cj) is an appropriate slope limiter, [Swe84]. There are many options 

for a limiter function. Some of the most usual choices are 

• MinMod (MM) limiter: 0(6*) = max(0, min(l, 6')), 

• VanLeer (VL) limiter: 0(^) = 

• Monotonized Central (MC) limiter: (f){9) = max(0, min((l + 6')/2, 2, 26*)), 

• Van Albada (VA) hmiter: ^{9) = f±g. 

The last three limiters have been shown to exhibit sharper resolution of discontinuities 
since they do not reduce the slope as severely as (MM) near a discontinuity. The TVD2 
reconstruction is second order accurate except at the local extrema where it reduces to the 
first order. A remedy is to consider the UN02 type reconstruction. 

The UN02 reconstruction is a linear interpolation which is second order accurate even 
at local extrema, [HOST]. The values W^, i , i are defined as 

= + = W^m - ^^m, (3.11) 



where 



2 



DiW = Wi+i - 2Wi + Wi_i, m{x,y) = ^(sign(x) + sign(?/)) min(|x|, \y\). 

Using either (TVD2) or (UN02) reconstructions the semidiscrete finite volume scheme 
(3.3) is formally second order accurate. 

In order to achieve higher order accuracy we also employ WENO type reconstructions 
for the values W^i, W^j^i- We implemented 3rd and 5th order accurate WENO methods 

(also referred to as WEN03 and WEN05, respectively) as they are described in [LOC94]. 
For the sake of simplicity we only present the WEN03 case. In order to compute the 
approximations and W^i, we first compute the 3rd order reconstructed values 
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We define the smoothness parameters 

and the parameters do = f , c^i = | and = di, di = d^, along with the weights 

OiQ tto ~ CKo ~ Oil 



where a, 



Wo 



ao + ai 



ao + "i 



ao + ai 



and e to be a small, positive number (in our computations we 



set e = 10^^^). Then the reconstructed values are given by the following formulas 



1 



r=0 



2 



r=0 



(3.12) 



3.1.4. Discretization of source terms. The finite volume discretization of the source term 
S{w) in (3.1) depends on the particular choice. On the other hand the resulting approxi- 
mation should preserve the upwind nature and the overall scheme should be well balanced. 
One possible discretization of the source term S{w) is given by: 



^ S{w) dx 



5 



i+4 



S 



i+i 



(3.13) 



3.1.5. Fully discrete schemes. Equation (3.3) is an initial value problem and can be dis- 
cretized by various methods. In our case we use a special class of Runge-Kutta methods 
which ensure the TVD property of the finite volume scheme, [S088, GSTOl, SR02]. 

Let At be the time step and let t""*"^ = + At, n > be discrete time levels. Assuming 
at t" the approximations {W^"}, i € Z are known then W[^~^^ are defined by 



n+l 



r. 



+ AtJ2b,Sr^^ 



+ At^ajeSi'^ 



(3.14) 



where J^^'\ 



The set of constants A 



(6i, . . . ,bs) define an s— stage Runge-Kutta method. The following tableau are examples 
of explicit TVD RK-methods which are of 2nd and 3rd order respectively 












1 





1 


1 


1 




2 


2 







1 
i i 

! t 2 



(3.15) 



In our computations we mainly use the three stage 3rd order method. 
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3.2. Finite volume schemes for dispersive models. To construct the finite volume 
schemes for the dispersive models the main idea is to rewrite the governing equations or 
systems in a conservative like form and discretize the resulting conservation laws using the 
aforementioned framework. One can use any of the numerical fiuxes, J-"*", J-''^^, J^^^ and 
reconstruction techniques TVD2, UN02 or WENO. Temporal discretization is based on 
the TVD-Runge-Kutta methods, (3.15). 

3.2.1. Boussinesq systems with flat bottom. Boussinesq systems (1.2) can be rewritten in 
a conservative like form as follows: 

(I-D)v,+ [F(v)L + [G(v)L = 0, (3.16) 

where v = (//, u)^, F(v) = ((1 + //)«, + ^m^)^, G(v) = {aUx^,cr]^xV , and D = 
diag{bd'^,ddl). The simplest discretization is based on the average fiuxes J-"™ for F and 
Q"^ for G. For the other two choices of the numerical fiux J-" the evaluation of a Jacobian 
is needed. Let A denotes the Jacobian of F, then 



A 



U 1 + T] 
1 U 



with eigenvalues Aj = m ± a/1 + r/, z = 1, 2. It is readily seen, since F is a hyperbolic fiux, 
that A can be decomposed as A = LKR thus for the characteristic fiux J^^^ we have with 



^ = Si = sign(Ai), z = 1, 2 



AiW, V) 



|(S1 + S2) + fJ'li^l - S2) 



For evaluating the numerical fiuxes J-", Q simple cell averages or higher order approxima- 
tions such as UN02 (3.11) or WENO (3.12) can be used. 

Remark 1. The discretization of the elliptic operator!) is based on the standard centered 
difference. This is a second order accurate approximation and it is compatible with the 
TVD2 and UN02 reconstructions. For higher order interpolation we modify the elliptic 
and flux discretization. Indeed, the finite volume scheme is modified as 



d 
di 



Vi_i + lOV, + V,+i Vi+i - 2Vi + V, 

{b,d)- 



12 ' ' ' Ax^ 

^ 12 ~ ' 

where Tii = ■^(J-'j+i — J^i-i) + -^(^j+i — resulting in a high order accurate ap- 

proximation. Thus in the WEN03 case a global third order accuracy is observed, while for 
WEN05 interpolation, we profit only locally by the 5th order accuracy of the reconstruction, 
cf. Section 4-T 

Remark 2. In the sequel for the discretization of the dispersive term G we use mainly the 
average numerical flux Q"^ defined as G^i = (a, c)|(Yj + Yj+i), where Yj = ■^^(Vj_|_i — 

2Vj + Vj_i). In case of higher order WENO reconstructions we use the average numerical 



12 



D. DUTYKH, TH. KATSAOUNIS, AND D. MITSOTAKIS 



flux based on the reconstructed values ofYi i.e. the flux G^^i = («-, c)|(Y^i + Y^^i), 
where i and i are reconstructed values of Y,- . 

3.2.2. Boussinesq system with variable bottom. We write system (2.15) in terms of depen- 
dent variables v := {H, Q)^ in the following conservative form 



[D(v,)] + [F(v)]. = S(v), (3.17) 



where 



-j-j^^ ^ I -^t j 18) 

We consider a uniform mesh and we denote by Hi, Ui and Di the corresponding cell 
averages. To discretize the dispersive terms in (3.18) we consider the following approxima- 
tions: 



1 
Ax 



X. 1 
'-2 



l + ^{Hxf~^HHx 



Q dx 



^ . 1., (3.20) 

Ax ^ 3 ^ 3 ' Ax2 ^ ^ 

The aforementioned discretizations lead to a linear system with tridiagonal matrix denoted 
by L that can be inverted efficiently by a variation of Gauss elimination for tridiagonal 
systems with computational complexity 0{d), d-being the dimension of the system. We 
note that on the dry cells the matrix becomes diagonal since Hi is zero on dry cells. For the 
time integration the explicit third-order TVD-RK method, (3.15) is used. In the numerical 
experiments we observed that the fully discrete scheme is stable and preserves the positivity 
of H during the runup under a mild restriction on the time step At. 

Therefore, the semidiscrete problem of (3.18) - (3.19) is written as a system of o.d.e's in 
the form 

L.v., + ^(J-._,i-J-._|) = ^Si, (3.22) 

where Lj is the i— th row of matrix L and J%- , i can be chosen as one of the numerical flux 
functions mentioned in the previous sections. In the sequel we will use the KT and the CP 
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numerical fluxes. In this case the Jacobian of F is given by the matrix 

1 
oH-{Q/Hf 2Q/H^ 



A 



and the eigenvalues are Ai^2 = Q/H ± y/gH. Therefore, the CF numerical flux takes the 
form 



i+4 



F(Vf J + F{Vf ,) F(Vf^ J - F(Vf 



(3.23) 



where /x = (/ii,/i2)^ are the Roe average values, 



Hi 



i+ 



i+4 



1^2 



i+4 



and 



S2(/J2+c)-Sl(/^2-c) SI-S2 
2c 2c 

(s2-gl)(A'i-C ) Si{^l2+c)-S2{lJ.2-c) 
2c 2c 



, c = ^/gJH, Si = sign(Ai 



(3.24) 



In order to guarantee the positivity of the reconstructed values -ffj+i on the cell interfaces 
we employ the well balanced hydrostatic reconstruction algorithm, [ABB"'"04]. Here we 
briefly recall the great lines of this reconstruction algorithm. 



In the cell Ci we compute flrst the reconstructions Vj^,. and Vi^i at {i + h) and {i 



respectively using either TVD2 or UN02 with MinMod limiter. Moreover, we compute in 
the same way the values rji^i and rji^r of the free surface elevation rji = Hi — Di. Now we can 
deduce the values Di i = Hi i — rji^i and Di r = Hi r — r]i^r- Letting D-^i = min(D, A,;) 
we compute 



H 



R 

i+l 



max(0, Hi^r + A,r - A+i)^ H^^i = max(0, + A+i,« - A+i 



and we deduce conservative reconstructed variables 



H,,,Ui,. 

1-1-2 



V 



H^.Ui+ii 



(3.25) 



(3.26) 



Then the term S, can be written as Sj = S^, i + i + Scj, where 



i+i 











and 



9 2 



Dir, 



Numerical experiments show that the resulting scheme is well-balanced even for Boussinesq 
system of equations. 
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3.2.3. Boundary conditions. In the case of Bona-Smith type systems with flat bottom we 
consider herein only the initial-periodic boundary value problem which is known to be 
well-posed [ADM09]. 

In case of the modified Peregrine's system with an uneven bottom we use reflective 
boundary conditions. We note that for the classical Boussinesq system posed in a bounded 
domain / = [&i,62]; one needs to impose boundary conditions only in one of the two 
dependent variables, cf. [FP05]. In the case of reflective boundary conditions it is sufficient 
to take u{bi,t) = u{b2,t) = cf. [AD]. In [AD] it was also observed that during solitary 
waves reflection the derivatives rjxihi.t) = rix{b2,t) — )■ 0, while for other wave types these 
derivatives remained very small. 

In our case we consider analogous reflective boundary conditions taking the cell averages 
of u on the first and the last cell to be Uq = Un+i = 0. We don't impose explicitly boundary 
conditions on H. The reconstructed values on the first and the last cell are computed using 
neighboring ghost cells and taking odd and even extrapolation for u and H respectively. 
These specific boundary conditions appeared to reflect incident waves on the boundaries 
while conserving the mass. 

4. Interactions of solitary waves 

For the Boussinesq system (1.2) we present initially results demonstrating the accuracy 
of the finite volume scheme. We study the propagation as well as the interaction of solitary 
waves. In particular we consider head-on and overtaking collisions. 

4.1. Accuracy test, validation. We consider the initial value problem with periodic 
boundary conditions for the Bona-Smith systems (2.7) with known solitary wave solutions 
(2.8) - (2.9) to study the accuracy of the finite volume method. We fix 6^ = 8/10 in the 
system and an analytic solitary wave of amphtude rjo = 1/2 is used as the exact solution 
in [—50, 50] computed up to T = 200. The error is measured with respect to discrete 
and scaled norms i?^, E"^, namely 

Elik) 

where = {Uj'}i denotes the solution of the fully-discrete scheme at the time t'^ = kAt. 
The expected theoretical order of convergence was confirmed for all finite volume meth- 
ods presented above. Three indicative cases, demonstrating the order of convergence, are 
reported in Table 1 : a) for the average flux, b) for the KT flux with TVD2 reconstruction 
using the minmod limiter and c) for CF flux with WEN03 reconstruction. The order of 
convergence for the WEN05 method cannot be obtained since a 4th order discretization 
is used for the elliptic operator. 



\u - f/^|U,oo/||f/°IU,,oo, \\u - f/^|U,oo = max \u{xi) - f/f | 
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(a) Average Flux 



(b) KT-TVD2(MinMod) 



Ax 


Rate(E2) 


Rate(E^) 




Ax 


Rate(E2) 




0.5 


1.910 


1.978 




0.5 


2.042 


2.032 


0.25 


1.910 


1.954 




0.25 


2.033 


2.029 


0.125 


1.923 


1.937 




0.125 


2.026 


2.023 


0.0625 


1.936 


1.941 




0.0625 


2.021 


2.019 


0.03125 


1.946 


1.948 




0.03125 


2.017 


2.016 



(c) CF-WEN03 



Ax Rate(E^) Rate(E^ 



0.5 


2.976 


2.975 


0.25 


3.017 


3.022 


0.125 


3.031 


3.044 


0.0625 


3.042 


3.059 


0.03125 


3.051 


3.073 



Table 1. Rates of convergence. 



We also check the preservation of the invariant (2.3) by computing its discrete counter- 
part: 



If = ^AxU,2 + [(l + r/,Kf-c 



Vi+i - Vi 
Ax 



1 2 



Ui+1 



Ui 



Ax 



(4.1) 



as well as the discrete mass Jq 



Ax^-rji. Figure 1 shows the evolution of the amplitude 
and the invariant of the solitary wave up to T = 200. The comparison of various 
methods is performed. We observe that the UN02 reconstruction is more accurate while 
KT and the CF schemes show comparable performance. (We note that the invariant /q = 
1.932183566158 conserved the digits shown for all numerical schemes. In this experiment 
we took Ax = 0.1, At = Ax/2.) Figure 1 (a) and (c) show the evolution of the amplitude 
of the analytical solitary wave of the Bona-Smith system (6'^ = 8/10) and of the solitary 
wave produced by the solution of the analogous ordinary differential equations system of 
the classical Boussinesq system respectively. In the case of the classical Boussinesq system 
we took Cs = 1.2 and we used the method described in [DM08]. 

The well balanced property of the finite volume schemes is also verified numerically. We 
consider a uniform shore, cf. Section 5, including a wet and dry region. The bottom is 
also modified by adding a small parabolic type hump located at x = 40. We tested the 
steady state preservation of all fluxes and possible reconstruction techniques. The results 
are similar in all cases. In Figure 2 we present the case of FC flux along with UN02. 



4.2. Head-on collisions. The head-on collision of two counter-propagating solitary waves 
is characterized by the change of the shape along with a small phase-shift of the waves 



16 



D. DUTYKH, TH. KATSAOUNIS, AND D. MITSOTAKIS 



O o ^ 



0.46- 
0.45 
0.44 



V 


CF-UN02 


□ 


CF-TVD2 





CF-WEN03 


o 


average 



O C) 



fl 1.2 



V CF-UN02 

□ CF-TVD2 

CF-WEN03 

o average 



20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200 

t 

Th 



(a) Evolution of rj amplitude (Bona-Smith) 

0.48 1 , , , 



(b) Evolution of /f (Bona-Smith) 



0.47 



0.46 
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0.41 
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V CF-UN02 

□ CF-TVD2 
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(c) Evolution of 77 amplitude (classical Boussinesq) 

Figure 1 . Preservation of the solitary wave amplitude and conservation of 
the invariant : G"" flux with Minmod limiter 



as a consequence of the nonlinearity and dispersion. These effects have been studied 
extensively before by numerical means using high order numerical methods such as finite 
differences, [BC98], spectral and finite element methods [ADMIO, DDLMM07, PDOl] and 
experimentally in [CGH+OG]. In Figure 3 we present the numerical solutions of the BBM- 
BBM system (2.6) and the Bona-Smith system (2.7) with 6^ = 9/11 (in dimensional and 
unsealed variables) along with the experimental data from [CGH+06]. The spatial variable 
is expressed in centimeters while the time in seconds. The solutions were obtained using 
the CF-scheme with UN02 and WEN03 reconstruction using Ax = 0.05 cm and At = 0.01 
s. For this experiment we constructed solitary waves for Boussinesq systems by solving 
the respective o.d.e's system in the spirit of [BDM07] such that they fit to experimentally 
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-20 20 40 60 80 100 20 40 60 80 100 

X X 

(a) t — 100 (b) magnification of tlie solution at t = 100 



Figure 2. Steady state preservation 

generated solitary waves before the collision. The speeds of the right and left-traveling 
solitary waves are Cr^s = 0.854 m/s and q^s = 0.752 m/s respectively. 

We observe that Boussinesq models converge to the same numerical solution with all nu- 
merical schemes we tested. A very good agreement with the experimental data is observed. 
The maximum height predicted by the numerical solution during the collision process is 
slightly higher in the case of the BBM-BBM system but the difference is negligible within 
the specific experimental scale. Furthermore, we observe similar underestimation of the 
maximum amplitude of the colliding waves compared to the experimental data, [CGH+06]. 
This discrepancy might be explained by a possible "splash" phenomenon during the colli- 
sion reported also earlier by T. Maxworthy, [Max76]. After the collision we observe that 
the phase shift of the solitary waves is the same in both numerical and experimental data, 
while the shape of the experimental solitary waves were not stabilized due to interactions 
with other small amplitude dispersive waves. We note that after the head-on collision of 
the waves small amplitude dispersive tails were developed, [BC98, ADM 10, DDLMM07]. 

The discrete mass for the Bona-Smith system is Jq^ = 0.0059904310418 and for the 
BBM-BMM system is Jq^ = 0.0059199389479 for all fluxes and reconstructions used. The 
variances in J{* are mainly due to different types of reconstruction and not to the choice of 
numerical fluxes. In Table 2 these values are reported. 

4.3. Overtaking collisions. The overtaking collision of two solitary waves similarly to 
the head-on collision incorporates nonlinear and dispersive effects. Overtaking collision 
has been studied recently in the case of bidirectional models in [ADMIO]. The interaction 
is similar to that of the unidirectional models but it was found that a new N-shape wavelet 
is generated during the interaction. This wavelet is of small amplitude and travels in 
the opposite direction to solitary waves and its shape depends on the Boussinesq system 
in use. Furthermore, as it was observed numerically and experimentally in [CGH+06], 
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(c) t = 19.00956s 
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100 



(b) t = 18.80067s 
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■ Expehmental data 






1 









Bona-Smith eWll 

- - BBM-BBM 
■ Experimental data 



300 -300 



(d) t = 19.15087s 



Bona-Smith 8 =9/11 

- - BBM-BBM 
■ Experimental data 



(e) t = 19.19388s 



(f) t = 19.32904s 



Figure 3. Head-on collision of two solitary waves: — : BBM-BBM, : 

Bona-Smith (^^ = g/n)^ experimental data of [CGH+06] 
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Figure 3. (Cont'd) Head-on collision of two solitary waves. — : BBM- 
BBM, : Bona-Smith {9^ = 9/11), •: experimental data of [CGH+06] 



(a) Bona-Smith 







m-flux 


0.000944236 


UN02 


0.00094423 


TVD2 


0.00094 


WEN03 


0.00094423 



(b) BBM-BBM 







m-flux 


0.00092793 


UN02 


0.00092793 


TVD2 


0.00092 


WEN03 


0.00092793 



Table 2. Preservation of the invariant 



the interaction of two solitary waves during an overtaking collision is characterized by a 
mass exchange and not by a simple superposition of the solitary pulses. These pulses 
remain separate retaining two different maxima contrary to unidirectional models where 
they merge into a single pulse momentarily. 

To study this interaction we solve numerically the Bona-Smith system (2.7) with 6"^ = 
9/11. Following the same process as before two solitary waves were generated numerically 
with speeds Ci^s = 1-2 and C2,s = 1.4. We solved the system using all fluxes using UN02 
and WEN03 reconstructions with discretization parameters Ax = 0.01, At = 0.005 up 
to T = 600. During simulations we were able to observe the generation and propagation 
of a small N-shape wavelet. In all computations the invariants were Jg = 4.6098804880, 
if = 5.116 conserving the digits shown for all methods. 

In Figure 4 we present the interaction of two solitary waves. Figure 5 shows a magnifi- 
cation on the generation of a small wavelet along with the generation of dispersive tails as 
an effect of the inelastic interaction of two waves. In Figure 6 we observe that the overtak- 
ing collision is accompanied by an exchange of mass between pulses while both peaks are 
permanently present. The situation is different for unidirectional models where two pulses 
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(c) t = 400 



(d) t = 600 



Figure 4. Overtaking collision of two solitary waves of the Bona-Smith 
system with 6*^ = 9/11. 

merge during a few time-steps to travel as a single pulse. Up to the graphic resolution 
we could not observe any difference in numerical solutions between UN02 and WEN03 
reconstructions. 

4.4. Small dispersion effect. In this section we study the small dispersion effects on 
solitary waves of the classical Boussinesq system. The motivation for this study is the 
lack of theory supporting the breaking phenomena in Boussinesq systems contrary to the 
KdV equation. For this reason we employ the Boussinesq system with a = 6 = c = 0, 
d = 10~^ and we take the solitary wave of the Boussinesq systems (2.5) as an initial 
condition. In Figure 7 we present numerical results obtained with CF-UN02 and CF- 
WEN03 schemes. In these experiments we take Ax = 0.001 and At = Ax/2. The 
invariant Iq is 1.629096452537 preserving the digits shown during all simulations. The 
invariant /{* is not preserved by this model since the coefficient b is not equal to d. The 
oscillations generated in the case of the WEN03 reconstruction were larger compared 
to those generated by the UN02 reconstruction. Moreover, a new W-shaped wavelet is 
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Figure 5. Generation of a wavelet during the overtaking collision of two 
solitary waves of the Bona-Smith system with (P' = 9/11. 

generated traveling to the left. This small wavelet finishes by producing a secondary 
breaking very similar to that of the initial solitary wave. 

5. BOUSSINESQ SYSTEM WITH VARIABLE BOTTOM: RUNUP OF LONG WAVES 

The shallow water equations are routinely used to predict a tsunami wave runup and, 
subsequently, constitute inundation maps for tsunami hazard areas. One of the main 
questions we address in this study is whether the inclusion of dispersive effects is beneficial 
for the description of the wave/shore interaction. In this section we perform a comparison 
of numerical solutions to Boussinesq equations (2.15), shallow water equations (1.1) (solved 
by the same numerical method) and experimental measurements made by CE. Synolakis 
[Syn87] and J. A. Zelt [Zel91]. In these experiments we consider a bottom of the form. 
Figure 8, 

^ f —xtanB, x<cotB, 
—D[x) = < , ^ r, 

^ ' \ -1, x> cot ,9, 
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Figure 6. Overtaking collision of two solitary waves of the Bona-Smith 
system with 6"^ = 9/11: mass exchange process. 



In all experiments over a fiat bottom, D{x) = Dq, we use an approximate solitary wave 
solution of the following form: 



r]o{x) = Assech^ (A(x - Xq)) , Uo{x) = -c. 



r]o{x} 



Do + r]o{x) 



X 



3A, 



4(1 + ^ 



V6{l + As) y/il + As)log{l + As)-As 
V3 + 2As 



(5.i: 



(5.2) 



where As denotes the amplitude, Cg is the correct speed of the solitary wave propagation 
for classical Boussinesq equations. 

The first three experiments we tested are described in [Syn87] and deal with the runup 
of solitary waves on a shore with a mild slope of 1 : 19.85. The first is a non-breaking 
solitary wave with dimensionless and scaled amplitude As/Dq = 0.0185, the second one is 
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Figure 7. The small dispersion effect onto classical Boussinesq equations 
solutions. 
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Figure 8. Sketch of the problem setup. 
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a nearly breaking solitary wave with As/ Dq = 0.04, while the third experimental setup is 
a breaking solitary wave with As/Dq = 0.28. 

System (2.15) has some advantages over other asymptotically equivalent models with 
variable bottom. Namely, it shows excellent stability properties even for nearly breaking 
waves on the shore. However, for the simulation of strong breaking events, it is beneficial 
to include friction or dissipative terms taking into account turbulence generation. 

We also considered two experiments from [Zel91] concerning the runup of solitary waves 
on a shore with steep slope 1 : 2.74. These experiments shed some light on the differences 
between dispersive and non- dispersive models. 

Finally we consider a non-uniform sloping shore that contains a small pond demon- 
strating the capability of the modified Peregrine's system to handle simultaneously and 
correctly dispersive effects in two basins with different mean sea levels. 

In the sequel t denotes the dimensionless time scaled by the quantity gjDQ. Further- 
more, we denote by R the height of the last dry cell at a specific time instance. In our 
computations a cell is considered as dry if the total water depth Hi inside is less than 
5 ■ 10~^^. The quantity R will also be referred to as runup. The maximum runup will be 
denoted by -Roo- In all experiments the discretization parameters were taken to be equal 
Ax = 0.05, At = Ax/lO, unless otherwise mentioned. Further, we compute in all cases the 
discrete mass /q and show the preserved digits. We use the KT and CF schemes combined 
with the TVD2 and UN02 reconstructions. The CF-scheme appeared to be less dissipative 
and we emphasize the results of this method. 

5.1. Runup of a solitary wave on a gradual slope /3 = 2.88° with Ag/ Dq = 0.0185. 

We consider first the simplest case — the runup of a non-breaking solitary wave. In this 
experiment we take an initial solitary wave with the amplitude Ag = 0.0185, Dq = 1 
and Xq = 19.85 in / = [—10,70] and a mildly sloping shore 1 : 19.85. This specific 
solitary wave does not break [Syn87] and the solution remains smooth during the runup 
and the rundown processes. In Figure 9 we show several profiles of numerical solutions to 
Boussinesq and shallow water equations along with the experimental data of [Syn87]. We 
observe that both models converge to the same solution. The runup as well as the rundown 
in this experiment is predicted very well. The runup value R for both models is almost the 
same. The maximum runup is R^o ~ 0.085 for the Boussinesq system, while for NSWE 
is Roo ~ 0.088. The experimental value reported in [Syn87] is equal to i?oo ~ 0.078. In 
Figure 10 the runup i? as a function of time is represented. The discrete mass is preserved 
/q = 60.3667671231 conserving the digits shown for both models. 

5.2. Runup of a solitary wave on a gradual slope f3 = 2.88° with As/Dq = 0.04. 
We consider the same sloping shore as before. We study the runup of a solitary wave with 
amplitude Ag = 0.04, placed initially at Xq = 19.85 in J = [—10,70]. The solitary wave 
does not break during the runup phase. Breaking occurs during the rundown process as in 
experimental observations [Syn87]. Results of the numerical simulations are presented in 
Figure 11. In Figure 12 the evolution of the runup value is shown. The maximum runup 
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Figure 9. Solitary wave runup oi 
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a sloping shore: Ag = 0.0185 case. 



for the Boussinesq system is Roo ~ 0.20 and Roo ~ 0.21 for shallow water system. The 
experimental value reported in [Syn87] is Roo ~ 0.156. 

In Figure 13 we perform a comparison with tide gauge data (free surface elevation 
measured in [Syn87]) collected at 32.1 m from the still shoreline position. We observe 
again a good agreement between the dispersive and nondispersive models. The discrete 
mass is preserved, Jg = 60.5210181987 conserving the digits shown. 

5.3. Runup of a solitary wave on a gradual slope (/3 = 2.88°) with Ag = 0.28. 
Finally we present the stiffest case of a solitary wave with amplitude As = 0.28, placed 
initially at Xq = 19.85 in J = [—10,60]. This specific initial condition is characterized 
by the wave breaking phenomenon before even reaching the shoreline. Strictly speaking, 
in this case Boussinesq model is not valid unless a wave breaking mechanism is consid- 
ered, cf. [Zel91]. In this case the approximate solitary wave given by formulas (5.1)-(5.2) 
with As = 0.28 it is outside the range of vahdity of the specific system. We proceed by 
constructing numerically a more accurate solitary wave following the cleaning procedure of 
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Figure 10. Runup value i? as a function of time: Ag = 0.0185 case. 



[BC98]. (Note that since the hnearized system (2.15) coinsides with the hnearized classi- 
cal Boussinesq system we conclude that there exist classical solitary wave solutions of at 
least small amplitude, cf. [DM08]). We constructed a solitary wave with Ag = 0.28 with 
Ax = 0.1 and Aa; = 0.05 while At = Ax/10. 

Specifically, we consider the interval [—200,200], and the initial condition ?7o(a;) = 
0.5 exp(—(x — 150)^/25), uq = (with reflective boundary conditions). The initial condition 
is resolved into two solitary waves traveling in opposite directions. We observed that the 
left-traveling solitary wave at t = 88 was separated enough from the rest of the solution. 
This solitary wave is isolated by keeping the numerical solution in the interval [—200, —151] 
and setting it equal to zero to the rest of the interval. The solution is then translated to 
the right at Xq = 19.85 and we let it propagate. We observe that the clean solitary wave 
propagates like the analytical solution of the Bona-Smith system with analogous behavior 
between the TVD2 and UN02 reconstruction. In Figure 14 we present the results for the 
KT scheme since the results of the CF scheme are completely analogous and no difference 
can be observed. In Figure 14 (a) the rightward traveling waves were reflected by the right 
boundary. 

In order to ensure the stability of the simulation and to study the runup, instead of 
smoothing, filtering or adding extra dissipative terms, we simply exclude the contribution 
of the term Qxxt in the vicinity of the shoreline (where Di < 0.3). Wave transition between 
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Figure 11. Solitary wave runup on a sloping shore: As = 0.04 case. 
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Figure 11. (Cont'd) Solitary wave runup on a sloping shore: Ag = 0.04 case. 
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Figure 12. Runup value i? as a function of time: Ag = 0.04 case. 



these two regions appeared to be smooth as one may witness on Figure 15. After this 
slight modification, the algorithm became more robust for large amplitude breaking waves 
without creating any unphysical oscillations. 

In this experiment friction appeared to play a significant role during the runup process, 
contrary to previous cases. The maximum runup computed without taking into account 
the friction of the bottom was far away from the experimentally measured values. For this 
reason, and only in this specific test case we included the empirical friction term (2.16) 
into the momentum conservation equation (2.15), with coefficient = 2 • 10"^. The 
friction term is discretized according to (3.13). This discretization preserves the positivity 
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Figure 13. Free surface elevation at x = 32.1 m. 

of all numerical schemes we tested. Mass conservation in this experiment was perfect 
Iq = 51.7504637472 preserving the digits shown. 

In Figure 15 we show the propagation of a breaking wave including its runup and run- 
down. We observe a significant difference between shallow water system and the dispersive 
model during the wave propagation. Discrepancies are present in the amplitude and in 
the phase speed simultaneously. However the dispersive model solution approximates bet- 
ter the measurements of J. A. Zelt [Zel91]. Nevertheless, we have to underline that the 
runup and rundown are fairly well described by both models. The maximum runup value 
according to the dispersive and nondispersive models is Roo ~ 0.47 which is in the range 
of [0.42,0.53] of the theoretical prediction of CE. Synolakis, [Syn87]. There is no single 
experimental value reported for the maximum runup in [Syn87] due to practical difficulties 
in generating a solitary wave of such large amplitude. Finally, we mention that the specific 
technique for handling the breaking wave leads to more accurate results for the rundown 
process than one presented in [Zel91]. 

5.4. Solitary wave runup on a steep slope (/3 = 20°). Now we present two experiments 
pointing out some further differences in solutions to dispersive and nondispersive models. 
These experiments were performed by J. A. Zelt, [Zel91]. We consider two waves in J = 
[—10,30] with amplitudes As = 0.12 and Ag = 0.2 initially located at Xq = 8.85 and 
Xo = 10.62 respectively. These waves propagate onto a steep sloping shore 1 : 2.74. A 
very fine grid of Ax = 0.01, At = Ax/100 is used to guarantee the accuracy and stability 
of simulations. 

As it was observed in [Zel91] both waves do not break during the runup but the second 
one generates a strong breaking event during the rundown. Friction does not play an 
important role in these experiments. Consequently, no friction term is included into the 
models. 
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Figure 14. Generation of a solitary wave with Ag = 0.28 applying cleaning 
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Figure 15. Solitary wave runup on a sloping shore: Ag = 0.28 case. 



Figure 16 shows the runup value i? as a function of time. We observe for both models that 
there is a phase lag compared to the experimental data. We believe that this discrepancy 
can be removed by changing the definition of the last dry cell. We also observe that shallow 
water equations over-predict the maximum runup and the minimum rundown while the 
Boussinesq model predicts correctly the extrema in both cases. 

Figure 17 shows the rundown of the solitary wave of amplitude Ag = 0.2 during breaking. 
One may observe that the experimental data consist of two curves due to the difficulty in 
measuring the surface elevation of the breaking wave due to 3D effects which become 
important. On Figure 18 the free surface elevation at a gauge is presented. The gauge 
is located 8.85 meters away from the still shoreline position. The reflected wave appears 
in both cases to be highly dispersive thus, the Boussinesq model provides a much better 
approximation. In the case Ag = 0.2, the mass remained equal to Jq = 29.770808175, while 
in the case case Ag = 0.12, Jq = 29.4861671693 conserving the digits shown. 
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Figure 15. (Cont'd) Solitary wave runup on a sloping shore: As = 0.28 case. 
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Figure 16. Runup value as a function of time. 
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Figure 17. Rundown of the wave with amphtude As = 0.2. 
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Figure 18. The amplitude at the wave gauge A 

5.5. Solitary wave runup on a gradual slope (/3 = 2.88°) with a pond. We repeat 
the experiment of Section 5.2 in / = [—10,50] with solitary wave amplitude equal to 
As = 0.05. However, we modify the bottom by adding a small pond over the shoreline 
described by the exponential function 0.1 e"^^"'"*^-' . The Boussinesq system preserves the 
correct dispersion characteristics for the waves reaching the pond. In Figure 19 we present 
the overall process. It is worth noting that after the pond was filled a breaking wave was 
reflected back. As the wave slides down, a small hydraulic jump appears. In the case of 
shallow water system this jump propagates as a shock wave due to hyperbolic character 
of equations. On the other hand, the Boussinesq system develops into an Airy type wave 
according to its dispersive characteristics. In Figure 20 we show the solution at two wave 
gauges located at x = —3.4 and x = 8 for both the dispersive and nondispersive models. 
The mass during the simulations is constantly equal to Jq = 40.5198087147. 

6. Conclusions 

In the present study we extend the finite volume framework, developed for hyperbolic 
conservations laws, to approximate solutions of dispersive wave equations. This type of 
equations arises naturally in many physical problems. In the water wave theory dispersive 
equations have been well known since the pioneering work of J. Boussinesq [Bou71] and 
Korteweg-de Vries [KdV95]. Currently, the so-called Boussinesq-type models become more 
and more popular as an operational model for coastal hydrodynamics and other fields of 
engineering. 

We extend the finite volume framework to dispersive models. We tested several choices 
of numerical fluxes, various reconstruction methods ranging from classical MUSCL type to 
modern approaches such as WENO. Various choices of limiters have been also tested out. 
Advantages of specific methods are discussed and some suggestions are outlined. 

For operational modeling of the wave runup we derived a new system which has some 
advantages over its classical counterpart. The new system together with proposed novel 
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Figure 19. Long wave runup on a shore with a pond. 
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Figure 19. (Cont'd): Long wave runup on a shore with a pond. 
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Figure 20. Evolution of the free surface elevation at two wave gauges. 
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discretization procedure are validated by extensive comparisons with experimental data of 
C.E. Synolakis [Syn87] and J.A. Zelt [Zel91]. 

We paid a special attention to the comparison of dispersive (Boussinesq) and nondis- 
persive (shallow water) models. Nowadays shallow water equations have become the 
model of choice for operational tsunami modeling including the inundation zone estima- 
tion [TG97, SBT+08]. The question of dispersive effects importance arises recurrently in 
the tsunami wave modeling community [KML05, DT07]. Our results show that shallow 
water equations are sufficient to predict maximum runup values. However, the dispersive 
effects can be beneficial for more accurate description of long wave propagation, runup and 
rundown. 
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